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MATRIX METHOD OF DETERMINING THE LONGITUDINAL-STABILITY COEFFICIENTS AND 
FREQUENCY RESPONSE OF AN AIRCRAFT FROM TRANSIENT FLIGHT DATA 1 

By James J. Donegan and Henry A. Pearson 


SUMMARY 

A matrix method is presented for determining the longitudinal- 
stability coefficients and frequency response of an aircraft 
from arbitrary maneuvers. The method is devised so that it 
can be applied to time-history measurements of combinations 
of such simple quantities as angle of attack , pitching velocity, 
load factor , elevator angle, and hinge moment to obtain the 
over-all coefficients. Although the method has been devised 
primarily for the evaluation of stability coefficients which are 
of primary interest in most aircraft loads and stability studies, 
it can be used also, with a simple additional computation, to 
determine the frequency-response characteristics. The entire 
procedure can be applied or extended to other problems which 
can be expressed by linear differential equations. 

INTRODUCTION 

The longitudinal characteristics of an aircraft are often 
related by a second-order linear differential equation in 
which the aircraft is assumed to have freedom in pitch and 
in vertical motion; changes in forward velocity are so small 
that they can be neglected. In the evaluation of tail loads, 
the coefficients of the differential equation and the elevator 
forcing function are generally assumed to be known and the 
response is to be determined. In the evaluation of gust 
problems the response and the coefficients are assumed to 
be known and the forcing function is to be determined. By 
analogy in stability and control work, it is desirable to 
determine the restoring-force and' damping-force coefficients 
from known forcing functions and responses. In case the 
damping is small enough to obtain the rate of decay (or 
logarithmic decrement) and period from the oscillation, the 
required damping and restoring coefficients are easily com- 
puted. Models employed in rocket-powered and drop tests 
can be and usually are so ballasted that such well-defined 
oscillations are obtained; however, the longitudinal oscilla- 
tions of piloted airplanes ordinarily are nearly critically 
damped and this analysis procedure cannot be applied. In 
any case, additional data and analysis are required to 
evaluate the control-effectiveness coefficients. 

Appreciable work has been done recently in the field of 
determining the frequency-response characteristics of air- 
craft in flight and evaluating the stability coefficients from 


the frequency-response data. In general, the methods for 
determining these relationships have been to impose actually 
prescribed motions such as unit steps, triangular pulses, or 
sinusoidal motions to the elevator by means of special equip- 
ment and then to measure the responses. The theoretical 
methods for reducing such data are usually tailored to fit 
the prescribed elevator motion. References 1 and 2 present 
methods of treating input and output data by Fourier 
analysis to determine the frequency response. Compared 
with the direct sine-wave input method of evaluating the 
frequency response, these methods require less special equip- 
ment and flight time at the expense of additional computa- 
tion. For the practical application of the Fourier transform 
method, certain restrictions are placed on the nature of the 
input and the resultant output motions: the motions must 
start from a trimmed steady-state condition and, at the end 
of the transient period, must approach either the original 
or the new steady or quasi-steady trim conditions. 

In view of the complications and limitations of existing 
methods of flight evaluation of stability coefficients and fre- 
quency response, development of a simple and less restricted 
flight test and associated analysis was considered desir- 
able. A matrix method for evaluating the longitudinal- 
stability coefficients of an aircraft directly from the input 
and output time histories corresponding to arbitrary con- 
trol motions has been derived in the present report. The 
frequency response and some of the stability derivatives 
may be evaluated once these coefficients are known. Al- 
though this method was derived to determine the second- 
order longitudinal response of an aircraft, it can be applied 
to other systems which can be approximated by second- 
order differential equations; extension of the method to 
higher-order linear systems is also possible. 

SYMBOLS 

Ai, A 2 combinations of aerodynamic parameters (see 
table I) 

b wing span, feet 

b t tail span, feet 

c chord, feet 

C h hinge-moment coefficient 



r Supersedes NACA TN 2370, '* Matrix Method of Determining the Longitudinal- Stability Coefficients and Frequency Response of an Aircraft From Transient Flight Data" by James J. 
Donegan and Henry A. Pearson, 1951. 
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C h( rate of change of hinge-moment coefficient with 

elevator angle (dCh/dS) 

C L lift coefficient ( LjqS ) 

C m pitching-moment coefficient of airplane without 

horizontal tail ( MbjqS 2 ) 

C mt pitching-moment coefficient of isolated horizontal 
tail surface 

g acceleration due to gravity, feet per second per 

second 

H hinge moment 

k v airplane radius of gyration about pitching axis, feet 

K empirical constant denoting ratio of damping mo- 

ment of complete airplane to damping moment 
produced by tail 
L lift, pounds 

m airplane mass, slugs (W/g) 

M pitching moment of airplane, foot-pounds 

n airplane load factor 

q dynamic pressure, pounds per square foot 



j S' wing area, square feet 

S t horizontal-tail area, square feet 

t time, seconds 

r true velocity, feet per second 

W airplane weight, pounds 

x, length from center of gravity of airplane to aero- 

dynamic center of tail (negative for conventional 
airplanes), feet 


KuK it -] 
Kz,K t , 
K S ,K S , 
K 1t K t , . 
K,J<C 


dimensional constants occurring in equations 
(see table I) 


a wing angle of attack, radians 

a, tail angle of attack, radians 

y flight-path angle, radians 

6 angle of pitch (a+y) 

5 elevator deflection, radians 

e downwash angle, radians (^- 

7] t tail efficiency factor (fff/?) 

<t> phase angle between incremental load factor and 

elevator deflection, degrees 
p mass density of air, slugs per cubic foot 

t dummy variable of integration 

« elevator angular velocity, radians per second 

The notations a and 9, a and 9, and so forth, denote single 
and double differentiations with respect to time. 
a bar over letter represents maximum value 

ja | bars on sides of symbol represent absolute value 



Figure l— Sign conventions employed. Positive directions shown. 


Matrix notation: 

|| || rectangular matrix 

[ ] square matrix 

{ } column matrix 

110111 integrating matrix (see table II) 

||A|| matrix defined by equation (24) 

|JA||' transpose of \\A | j 

Subscripts: 

i denotes row elements in matrix 

j denotes column elements in matrix 

t tail 


LONGITUDINAL EQUATIONS OF MOTION 


ELEVATOR MOTION 


In tins section the usual longitudinal equations of motion 
following an elevator motion are derived in such a manner as 
to obtain expressions between some of the simple combina- 
tions of variables which are measurable in flight: namely, 
angle of attack and elevator angle, pitching angular velocity 
and elevator angle, or load factor and elevator angle. The 
usual assumptions of linearity, small angles, no loss in air 
speed during the maneuver, and no flexibility are implied. 

As in reference 3, the differential equations of motion of 
an airplane due to a given elevator deflection may be written 
as (see fig. 1 for definitions): 


dC m 

da 


mjV-~- L - AaqS- 

da 


A aQ 


S 2 , dC L 

b + da t 


Aa^l — 


dCi 

dS 

de 

da 


Tj t qS t A5=0 


. x, de .x, 
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( 2 ) 
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By use of the definitions 

A9= A 7 -j-Aa ^ 
6=7- 
6 =7 + a 


( 3 ) 


equations (1) and (2) are reducible to the following second- 
order differential equation giving the relation between angle 
of attack and elevator angle: 

cc -}- A .2 Aa=A.3 A5 -|--Za4 5 (4) 


where the K’s are the constants for a given set of conditions 
and are defined in table I. The coefficient K represents an 
effective aerodynamic-damping coefficient; K 2 represents an 
effective aerodynamic restoring-force coefficient; K s and Ki 
represent effective elevator-control power coefficients. 

An alternative form of equation (4), expressing the relation 
between angle of pitch and elevator angle, may be obtained 
by inserting relations of equation (3) into equation (4) and 
noting from equation (1) that 


7 = Ai Aa-p AjA 5 

7 = A,ci -f A-2 5 

A7 = A, | *Aa dt+Atj A Sdt 


(5) 


where A\ and A 2 are combinations of aerodynamic param- 
eters defined in table I. The equation obtained after 
these preceding substitutions are made is 

e+KiB+KtM^W+Ktj'ASdt (6) 


where (see table I) 


K^Kt+KtAt+BUAi 

^AtKt+AiK, 


From the following definition for load-factor increment 


, r . dc L 

A n = — 7 — 

9 


it follows that 


a V 

A a-\ At AS 

da W S g 


Act- 


W/S , A, 

■wr 

da ® 

W/S . At . 

-dcr n ~A l 5 


W/S 

' d.C L n A, 

doc q 


A* •• 


( 7 ) 


(8) 


TABLE I.— DEFINITION OF CONSTANTS OCCURRING 
IN EQUATIONS 


Con- 

stant 

Definition 

Ai 

P V rdC Lt S t x? ( K de\ 

2 m[_da, k d V ‘ \^—'dcc) 


a 2 

pV 2 fdCi, S 2 , dC L , S t x t 

f (, de\ dC L K p&H) 

2m ( da k s -b ' da t >ft fed 

Lv da ) da v Tt 2 m Ji 

a* 

pV^rdC-L, S ffl , dC mi st dC Ll dC Lt KntpxtSt 1 

2m\_ d& nt V " r ds Vt b l k * da, dS 2 mkj J 

Ki 

d( K, S, 
dS Vi mV q 

-d-l 

dC L „ 

qS 

mV 

■d-2 

d ° L ‘ 

ds v ‘ qS ‘ 
mV 

K s 

Kz'T KiA.2 “f" K±Ai 


A 2 A 3 + AjAj 

Ai 

dC L 

da q K A K 

IF/S Ki ' g AlKi 

Ai 

dC_L 

IF/S Ki+ g A '- Kl 

K, 

E Ai 

g 

Ad 

Aj Z>C kx t ( de , l\ 
1 C ki ba t V\da' y !~) 

Ad 

k ' El Eljif i — d i_ d Ek l±\ 

2 _r C K, bat \ da da 2m yT(/ 

Ad 


Ad 

dC L 

da 9 Ad 
IF/S 3 


A ■* 

The — p 5 term in equation (8) was found to be small and 

Ai 

is omitted in the subsequent derivation. 

Substituting the results from equation (8) into equation 
(4) yields another form expressing the relation between 
measured load-factor increment and elevator angle as 


h + Kin T K-i An =K 7 ASAK s 8 


( 9 ) 
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where (see table I) K 7 and K s are now different forms of the 
effective control power coefficients. 


HINGE MOMENT 

The coefficients K\ to K% occurring in equations (4), ( 6 ), 
and (9) are those associated with the measured elevator- 
motion case. The use of the relation 

c ’‘-^ a5 +si> (10) 

gives the solution for A 8 as 

A5= ^bC h (° h ~b^ Aat ) (1 1} 

dS 


The increment in tail angle of attack to be substituted in 
equation (11) is given by 



so that 


dC L p S x t \ . x, /dt . 1 I 

da 2 m 1 \da y^7/ J 

( 12 ) 



d 


£i 

r 



dt dC L p S x, \ 

da da 2 m 


(13) 


In order to shorten the subsequent derivation for the hinge- 
moment case, the term K 4 S in equation (4) and its counter- 
parts in equations ( 6 ) and (9) are omitted. This effect, is 
usually small; however, each individual case should be exam- 
ined to sec whether the term warrants dropping. 

A substitution of the value of AS given by equation (13) 
into equations (4), ( 6 ), and (9) gives the following three 
differential equations for the same combination of variables 
with C h and its integral replacing AS: 

a+K 1 °a+K t 0 Aa=K i °C ll (14) 

e+K l i 'd+K 2 °Ae=K i 0 C h + 9 ~^° f‘c h dt (15) 

n+K 1 0 n+K 3 °An=K i °C k (16) 


INTEGRAL FORM OF EQUATIONS 

Although equations (4), ( 6 ), (9) and (14), (15), (16) could 
be used to evaluate the effective K coefficients from flight 
measurements of A a, and An together with measurements 
of elevator angle, stick force, or hinge moment, it is seen that 
several differentiations of the measured data would be re- 
quired. Inasmuch as a numerical differentiation process is 


inherent!}- more inaccurate than the. corresponding integra- 
tion process, the preceding equations are changed and 
rearranged so that either Aa, 9, or An, which are to be the 
measured values, appear as separate quantities on one side 
of the equation and the operations on these quantities 
appear on the other side. In integral form the rearranged 
equations are 



r Aa 
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f ASdrdt — 
Jo 




Addt — —A 

>0 

a 
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K6 
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An 
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dt+K, ° C f. 

A a drdt — K 3 ° 

P r a dr dt = 

— A a 






JO , / 0 

(20) 

K^e+Kz 0 

Cedt-Kt 

Jo 

C Chdt ~ 

. k 0 — 

-A 6 y 

f‘£ Kdrdt- 








(21) 

K ‘i 

"*{ 

An < 
0 

* + wf,T 

An drdt- 


j' ( \drdt= — 

An 


( 22 ) 


In principle to solve any one of these equations for the K 
coefficients, it is 011 I 3 - necessary to tabulate the recorded 
values of the two basic variables (for example, in equation 
(19) the values of An and AS) at a number of points t u t 2 , t 3 , 
and so forth along a given time history and to perform the 
indicated integrations from t =0 up to the time of the recorded 
value ti. A number of simultaneous equations containing 
the unknown K’s result which are then solved. The number 
of equations can vary from a minimum, in which the number 
of ordinates is equal to the number of unknown K's, to the 
case where there are more equations than unknowns. When 
the number of ordinates equals the number of unknown K’s, 
the usual methods of solving simultaneous equations may be 
used to obtain the IT’s; however, when there are more equa- 
tions than unknowns, a least-squares method is required to 
reduce the equations. Since the best average value of the 
li’s is obtained when many points along the time history are 
used, a least-squares procedure is generally preferable. 

Although the integration indicated in equations (17) to 
( 22 ) can actually be performed graphically from the time 
histories, it is deemed better to express the equations in 
matrix form in order to enable a complete numerical solution 
to be made. 
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MATRIX FORM OF EQUATIONS 

Since the derivation in matrix form for any one of equations (17) to (22) is the same as for any other equation, only 
equation (19). involving measured load factor and elevator angles, is used. In matrix form the system of simultaneous 
equations obtained from reading the time history of the load factor n against elevator angle § in an arbitrary pull-up may be 
written 


P An dt P V An irdt- pASdrdf- pASdi 

Jo Jo Jo Jo Jo Jo 


f "N 


r *\ 

— An L 
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-An, 

1 , , , , , . 


h. ^ 


V. * ^ 


(23) 


In shorter form this expression may be rewritten as 


r K^ 


Mil 


K, 

K-, 

IK, 


h = {-An,} 


(24) 


where the matrix ||A|| is in general a rectangular matrix; 
that is, for every time ft, one equation or one row of the 
matrix |[A|| is obtained. The individual elements of matrix 
||A|| are evaluated from the known values of incremental 
load factor and incremental elevator angle. As mentioned 
previously, the integration may be performed graphically 
but in the present case, use is made of the integrating 
matrices derived in reference 4. Thus, anv element in the 

rt 

rectangular matrix (equation (23)) such as Andt or 

Jo 

I I An d t dt may be expressed in matrix form as follows: 

Jo Jo 

f P- 


'An (ft j = ||CM| {An,} 


(25) 


f‘ (^Andrdt^liailjJjAn dtj 


The integrating matrix 11(7111 as derived in reference 4 is 
given in table II, with a time interval Af=0.1 second. It 
should be noted that a sufficient number of time intervals 
within the natural period being computed must be chosen to 
give a solution; usually the shorter the time interval chosen 
for the integrating matrix, the more accurate will be the 
final solution. 


After the elements of the matrix A (equations (23) and 
(24)) have been determined either by applying the inte- 
grating matrix or by graphical integration, the method of 
least squares is applied to the solution of the system of 
simultaneous equations. In matrix notation the least- 
squares solution involves multiplication of matrix A by its 
transpose A' so that equation (24) becomes 


[A'A] {£,} = { -A' An,} 


(26) 


where the matrix [A' A] would be a 4 by 4 matrix for equa- 
tions (18) and (19). Equation (26) can now be arranged to 
be solved directly for the iTs by multiplying by the inverse 
matrix [A'A] -1 so that finally 


r k i 

K, 
K 7 
K 4 


U [A'A]- 1 {—A' An,} 


(27) 


Alternately the system of simultaneous equations repre- 
sented by equation (26) can be solved for the values of K by 
any of the well-known methods of solving sets of simul- 
taneous equations, that is, by eliminating the variables or by 
using Crout’s method (reference 5). The derivation in 
matrix form of any of the other equations from (17) to (22) 
is similar to the plan given for equation (19) and, therefore, 
is not given. 

FREQUENCY RESPONSE 

As first derived by Cornell xAeronautical Laboratory 
(reference 6), the frequency response was measured by 
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TABLE II.— INTEGRATING MATRIX [C,] 


(Based on 0,1-see intervals] 


t 

0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

0.9 




0 

0 

0 

0 

0 

0 

o 

0 

0 

0 




- :■-? 

.1 

.041667 

. 0 G 6667 

-.008333 

0 

0 

0 

0 

0 

0 





.2 

.033333 

. 133333 

.033333 

0 

0 

0 

0 

0 

0 





.3 

. 033333 

. 133333. 

.075000 

.066667 

-.008333 

0 

0 

0 

0 





.4 

. 033333 

.133333 

. 066667 

. 133333 

.033333 

0 

0 

0 

0 





.5 

. 033333 

. 133333 

. 066667 

. 133333 

. 075000 

.066667 

008333 

0 

0 





.6 

. 033333 

. 133333 

. 066667 

. 133333 

.066667 

. 133333 

. 033333 

0 

0 





.7 

.033333 

.133333 . 

.066667 

. 133333 

. 066667 

. 133333 

.075000 

. 066667 

-. 00 S 333 




■ s - 

.8 

. 033333 

. 133333 

.066667 

. 133333 

.066667 

. 133333 

. 066667 

. 133333 

.033333 





.9 

. 033333 

. 133333 

.066667 

. 133333 

.066667 

.133333 

. 066667 

. 133333 

.075000 





1.0 

. 033333 

.133333 

. 066667 

. 133333 

.066667 

. 133333 

.066667 

. 133333 

.066667 





1.1 

.033333 

. 133333 

.066667 

. 133333 

.066667 

. 133333 

.066667 

. 133333 

.066667 





1.2 

.033333 

.133333 

.066667 

. 133333 

. 066667 

. 133333 

. 066667 

. 133333 

.060667 





1.3 

. 033333 

.133333 

.066667 

. 133333 

.066667 

. 133333 

.006667 

. 133333 

.066667 





1.4 

.033333 

. 133333 

.066667 

.133333 

.066667 

, 133333 

. 066667 

. 133333 

.066667 





1.5 

.033333 

. 133333 

.066667 

. 133333 

. 066667 

. 133333 

. 0 G 6 G 67 

. 133333 

.066667 





1.6 






+ 
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actually subjecting the airplane to sinusoidal elevator mo- 
tions of various frequencies by means of specially constructed 
apparatus. From these results the coefficients K j, K 2 , and 
so forth, which are significant in control and loads work, 
could be determined provided the equation of motion was 
assumed. 

In the present instance since the coefficients Ki and K 2 
are determined directly from the equation of motion, the 
corresponding relations are given so that the frequency 
response, which is significant in the design of stable autopilot 
systems, can also be determined. 

When a sinusoidal elevator motion has been assumed, 
then equation (9), omitting the minor effects of K$5, becomes 

h+Kin+K 2 An=K~o sin cot (28) 

where An is the load-factor increment and co is the angular 
velocity of the elevator. Since equation (28) is a linear 
equation with constant coefficients, the steady-state solutions 
are of the form 

n=n sin (cof+$) \ 

n—no> cos («f + <£) > (29) 

n ——nu> i sin (wf+#)/ 

By a substitution of these relations into equation (28) the 
following equation is obtained: 

— Tied 2 sin (wf+$) -bianco cos (&>£+<£) + 

K 2 n sin (o>f + <£)=iir 7 5 sin c d (30) 

which may be rewritten as 

n(K 2 — co 2 ) sin (cof+^+FTinco cos sin co t (31) 

or 

B sin (ut-{-<j>-\-t)=K 7 '5 sin cot (32) 

where 

B=K 7 8=n^(K 2 - co 2 ) 2 -}- (. K» 2 (33) 

and 

— f = ^=tan _1 (34) 

K-2 — CO" 

From equation (33) the amplitude ratio of load factor to 


elevator angle is seen to be 

|S|= Kl - (35) 

and the phase angle at various frequencies is given by 
equation (34). 

In the present case the values of K h K 2 , and Iv r would 
have been derived from the flight measurements and the 
values of co would be assigned. 

For the measured hinge-moments case the values of Kf 
and K 2 would be used instead of K\ and K 2 , and so forth. 
The complete frequency-response relations and transfer 
functions including all derivatives and integral of 5 for 
equations (4), (6), and (9) are given in the appendix. 


DETERMINATION OF AERODYNAMIC DERIVATIVES 


The various K coefficients determined from the meas- 
ured values may be termed effective coefficients and include, 
to some extent, effects of some nonlinearities, elasticity 
and effects of other variables which are omitted in the 
usual analysis. In addition, as may be seen from table 
I, the fTcoefficients are combinations of various quantities 
involving known geometric qualities, the conditions of the 
problem as well as aerodjmamic derivatives. The stability 
coefficients given in table I are expressed in a form suitable 
to loads work. In usual stability calculations, these co- 
efficients are generally expressed in a simpler form where 
the number of aerodynamic variables are reduced and, as 
a result, the coefficients are more easily approximated. 

dC L dC m de dC L , 

A total of 10 aerodynamic variables 


fa 


da ’ da 


dC L , dC, 


da t 


dS 


da da AT r , 

-jt-i rr - — , , v t , and K appear m the definitions of 

d5 da, OS 


the coefficients of table I. Although all the aerodynamic 
derivatives cannot be determined directly from the four basic 
coefficients (namely, K l; K 2 , Jv 3 , and K t ), engineering approxi- 
mations of the more significant derivatives can be obtained 
if values are assigned to either some of the more accurately 
known derivatives or to those factors having least influence 
on the problem. 
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Tlie factors having the least influence on the problem are 


K, Vt: 


and the derivative 



which, respectively, allow 


for the contribution of wing-fuselage damping, tail efficiency, 
and moment due to tail camber to which average values 
can be assigned. A representative value of K is 1.25. 
Representative values of tj, range from about 1.2 to 0.8 
with the higher limit applying to propeller-driven airplanes 
operating at low speed and full power and the lower limit 
applying at high speed with the propeller braking. An 
average value for jets or at zero thrust for propeller-driven 

dC m . 

airplanes is about 0.9. A representative value of , - 


can be obtained from existing wind-tunnel data or by using 
theoretical methods; —0.5 is an average value for tail 
surfaces. 

Since, as may be seen from table I, K t is directly propor- 


tional to 


<15 ’ 


an effective value of this derivative can 


be determined directly from the definition of fu. 

In order to determine consistent values of the remaining 


significant aerodvnamic derivatives 

da 


dC_ m 

da 


de d ( 'l, 
da da. 


dC\ 

da, ’ 


and further values must be assigned to several of the 

remaining derivatives. The derivatives chosen would natu- 
rally be those for which values could be obtained from other 
sources with the greatest degree of accuracy. 


EXAMPLES 

In order to illustrate the foregoing method as well as the 
consistency of results obtained with different sets of instru- 
mentation, typical examples are given using data obtained 
from three flights (referred to as flight 1, flight 2, and flight 3) 
of a high-speed medium jet bomber. For flight 1 , the method 
of a computation is obtained in sufficient detail to enable a 
reader not too familiar with the mathematical details to 
reproduce similar results. Flight 1 is further divided into 
case I where data for An and A 5 are used and case II where 
data for 9 and A5 are used. References 7 and 8 may be 
consulted for introductory discussions of least-squares and 
matrix methods. 

Figure 2 shows the measured time histories of velocity, 
altitude, incremental elevator displacement, incremental load 
factor, and incremental pitching velocity obtained during a 
push-down pull-up maneuver. By means of the values from 
figure 2, increments in load factor and elevator angle at 
0.1-second intervals have been tabulated in columns 2 and 3 
of table III. The. elements of the A matrix (equations (23) 
and (24)1 are given in columns 4 to 7 of table III. Each 
element in these columns has been determined by performing 
the indicated integrations on the results given in columns 
2 and 3. In this instance the integrations have been 
performed by use of the previously mentioned integrating 
matrix derived in reference 4. This method is particularly 



FKsURE 2.— Time histories of velocity, altitude, incremeatal elevator displacement, incre- 
mental load factor (computed and measured), and incremental pitching velocity (computed 
and measured) for Sight 1 at Mach number O.iO. 


suitable when automatic computing machines are available. 

The elements of matrix A (equation (23)) which are given 
in columns 4 to 7 of table III indicate that with the A t 
spacing used there are 23 equations involving the four 
unknown values of K. In order to obtain the least-squares 
solution of these equations, the transpose [[A|j' of matrix |[A|[ 
is required. The transpose matrix is obtained by inter- 
changing the rows and columns of matrix |[A|j. 

The product of the 4-row, 23-column transpose matrix by 
the 23-row, 4-cohunn original matrix yields the 4-row, 
4-column matrix in the coefficients of K t . The resulting four 
simultaneous equations are then solved by any of the well- 
known methods of solving sets of simultaneous equations. 

By performing the preceding operations, the following 
values of IF were obtained from the data listed in table III: 

K, K 2 K- K$ 

3.314221 7.339706 —119.553905 5.819025 

In order to show how well these computed values of K 
represent the original data, they have been reinserted into 
equation (19) along with the measured values of A5 to deter- 
mine calculated values of An. The computed curve is given 
by the dashed line in figure 2 of the plot of An against t. 
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TABLE III.— TABULATED VALUES FOR FLIGHT 1 


. ... . 

1 

2 

3 

4 

•5 

6 

7 

Time, t 
(sec) 

Acceleration 

increment, 

An 

Elevator angle 
increment, AS 
(radians) 

X And/ 

f f A ndrdt 

JO Jo 

— f f A& dr dt 

Jo Jo 

. 

JO ASdt 

0 

0 

0 

0 

o 

0 

0 

.1 

.054 

.046687 

.004050 

.000225 

-.000091 

-. 002465 

.2 

-.054 

. 077666 

.005400 

.000720 

-. 000622 

-.OOS8I4 

.3 

—.111 

.082902 

-.002133 

. 00096V 

001903 

-.016868 

.4 

-.254 

.085084 

-.019667 

-.000040 

004008 

-.025592 

. 5 

-.444 

.089010 

-.054950 

003632 

-.006969 

— . 033990 

.6 

-.588 

.093810 

106933 

-.011587 

010821 

-.043124 


-.784 

.098610 

— . 175658 

-. 025559 

-. 015610 

-. 052741 

.8 

-.965 

. 103846 

263233 

-.047347 

-.021386 

-.062860 

.9 

-1.122 

. 106900 

-.367483 

-.078747 

-.028197 

— . 073405 

1.0 

-1.291 

. 109081 

-.488033 

-. 121386 

036076 

-.084211 

1.1 

-1. 462 

. 108645 

-.626166 

— . 176978 

-.045041 

-. 095094 

1.2 

-1.575 

. 108645 

— . 778500 

-. 247093 

-.055094 

-.105954 

1.3 

-1.704 

. 108645 

943233 

-.333111 

-.066233 

-. 116830 

1.4 

-1. 739 

. 107336 

- 1.116166 

436013 

-.078457 

-.127610 

1.5 

-1.837 

.053668 

-1. 296083 

-. 556599 

-.091666 

-. 135C61 

1.6 

-1.8Q1 

.003491 

- 1 . 479099 . 

-. 695333 

-.105416 

13S490 

1.7 

-1. 700 

. 004800 

-1. 654399 

-.852104 

-.119294 

138995 

1.8 

-1. 569 

-.004799 

—1.818099 

-1. 025826 

133202 

-. 1390S6 

1.9 

-1. 305 

-.017017 

-1.961174 

-1. 214978 

-.147063 

-.137970 

2.0 

- 1.116 

-.026179 

-2. 081599 

-1.417305 

-.160760 

135785 

2.1 

-.869 

-.036651 

-2. 181332 

-1. 630682 

-.174186 

132600 

2.2 

-.564 

-.041887 

-2. 253466 

-1.852652 

-. 187254 

-. 128629 

2.3 

-.315 

-.045814 

-2. 296799 

-2.0S034! 

-. 199905 

-.124284 


The same process as was used for the relat ions of An and AS 
was also applied to the relations of 6 and AS shown in figure 2. 
The tabular material corresponding to table III is not 
included; the values of K obtained, however, were as follows: 

K\ I\. 2 Hs He 

3.13167 8.4123 -7.6212 -12.1967 

These values of K when reinserted into equation (18) 
resulted in the computed curve of d given by the dashed 
curve of figure 2. 

In addition to the preceding computations, several push- 
down pull-up maneuvers, made under similar conditions of 
altitude, weight, and center-of-gravity positions, were ana- 
lyzed to obtain the variation of several of the computed K ’ s 
with Mach number. In this analysis only, the measurements 
of An and AS were used. The results obtained for three Mach 
numbers are shown in figure 3. The short parts of the curves 
shown are the expected variations in the K' s. Table I shows 
that Iv, should vary linearly with speed and the other values 
of K should vary parabolicallv. . The curves shown are 
merely guides adjusted to pass through zero and through the 
value of K at the 0.45 Mach number point. 

The values of K { and K 2 shown in figure 3 were also inserted 
into equations (34) and (35) to determine the corresponding 
curves of frequency response. The results are given in 
figure 4. 

In addition the values of K u K 2 , and frequency response for 
case I have been computed by using the definitions of table I 
and aerodynamic derivatives obtained from wind-tunnel 
tests. These results are also shown in figures 3 and 4. The 
aerodynamic derivatives were listed in an unpublished report 
by the North American Aviation, Inc. and were obtained in 
the Southern California Cooperative Wind Tunnel. 


DISCUSSION _ _ 

If only the frequency response is desired, it can be deter- 
mined without recourse to the equations of motion; however, 
if the stability coefficients are desired, it will be necessary to 
use the equations of motion as has been done in the present 
report. For either case several mathematical methods are 
available (references 1, 2, and 6) to obtain these required 
quantities and all methods, if carried far enough, should 
yield similar results. Thus, the present method is basically 
no more accurate than any other method ; however, it has the 
advantage of simple instrumentation and experimental 
procedures but may require more extensive computation. 

As with other methods where linearity is a basic assump- 
tion, most consistent results are to be expected when the 
maneuvers are confined to the angle-of-attack region where 
linearity exists. In order for the outlined mathematical 
procedures to succeed, the maneuvers should cover as much 
of the linear range as possible in a short period of time and 
the portion of the maneuver considered should be confined 
to that portion where the integrals are increasing. This 
practice insures that the elements of the original matrix A 
are all different and that the subsequent least-squares 
matrix [A' A] is not ill-behaved. Enough of the response 
time history should be taken to cover a good portion of the 
natural period of the system. A point worth noting in 
connection with the use of the equations is that zero time is 
assumed as being at the start of the maneuver when the 
airplane is in steady flight. Since the present method is not 
restricted by the final condition, it offers the possibility of 
performing an analysis on fragments of curves with the result 
that any variations in the constants may be determined. 
In such an analysis two possibilities occur: (l) where the 
fragments considered start from a fixed initial condition and 
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.4 

Mach numcer 

Fl'R'RE 3.— Variations of K\ and Ki with Mach number. 


become successively longer, and (2) where the fragments are 
taken as consecutive. In the first case, the present method 
may be applied without any modification; in the second case, 
the equations must be altered to introduce the initial condi- 
tions for each fragment. These possibilities have not, 
however, been explored. 

In the derivation given herein, lag in downwash has been 
included (see equation (2)) but unsteady lift effects have not. 
References 9 and 10 show that for the present purposes the 
inaccuracy of omitting unsteady flow effects, except down- 
wash lag, is probably no greater than the inaccuracies in the 
original assumptions or of the experimental data. 

Other terms and other combinations of measurements 
might have been included in the derivations given- — for 
instance, the equations are readily adapted to measurement 
of tail load and either airplane load factor, airplane angle of 
attack, or pitching angular velocity. Additional terms may 
have been included to account for flexibility. Also it is 
possible, as for example in the case of the hinge-moment 
relations, to include additional terms to account for elevator 
moment-of-inertia effect, rate of elevator motion, and so 
forth in order to make the methods more inclusive. The 
inclusion of these further terms, however; generally requires 
additional Tv’s to be evaluated and would only be justified 
when the assumptions implied in the basic equations of 
motion can be more closely approached and when the 
accuracy of measurements is high. Although the method 
had been applied herein to second-order differential equations, 



a, radians) sec 

Figure 4.— Airplane frequency response. 


it may be extended to higher-order equations with the limita- 
tion that too many integrations destroy the conditioning of 
the equations used in determining the coefficients (equation 
(26)) and make the equations difficult to work with. 

The results of the sample computations in which two 
different sets of instrumentation were used indicate an 
average difference between the respective K coefficients of 
about 10 percent. The use of a least-squares method 
permits calculation of a probable error, which is an indication 
of how well the second-order system and the coefficients 
(computed on the basis of 0.1-second time intervals) fit the 
data. The expression used in computing the probable error 
is 

# ( = 0.6745^/^V^ 

where Bu is the main diagonal term of [AT4] -1 , E is the 
difference between the computed and measured value of the 
variable, N is the number of cases considered in the least- 
squares procedure, and k is the number of variables deter- 
mined. This probable error has been calculated for case 
I and case II and indicates an error of +0. 3 in Ki and ±0.5 
in Ki for the computations in which the accelerometer 
measurements were used. These values are contrasted with 
probable errors of ±0.1 and ±0.3 for the pitching-angular- 
velocity measurements. These probable errors are asso- 
ciated with the very small differences between the solid-line 
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and dashed-line curves shown in figure 2 . Greater accuracy 
may be obtained by increasing instrument accuracy, record- 
reading accuracy, and correcting original data for instrument 
errors. Further accuracy in the method may always be 
attained by using smaller time intervals. 

The results shown in figure 3 for the three flights investi- 
gated give some idea of the scatter to be expected between 
runs as well as the variation of the coefficients K-, and K> with 
Mach number. As might be expected from the definition, 
A, is seen to vary linearly with Mach number with little scat- 
ter. On the other hand, the values of A 2 either indicate a 
linear variation with Mach number or a scatter about the 
expected parabolic variation. 

The computed values of K x and A 2 (fig- 3) obtained from 
the wind-tunnel data are in fair agreement with the flight- 
test values. For many engineering purposes this agreement 
may be adequate and probably typical of what might be ex- 
pected if wind-tunnel data were used at the design stage. 
Since all of the A’s are defined in table I, the dynamic longi- 
tudinal characteristics of an aircraft may be estimated in the 


design stage by computing the A’s and inserting the values 
in the frequency-response relations given in the appendix. 

CONCLUDING REMARKS 

A matrix method has bfcen presented for determining the 
longitudinal-stability coefficients and frequency response of 
an aircraft from an analysis of arbitrary maneuvers in which 
simple instrumentation is used. Errors in instrument accu- 
racy and probable errors due to the use of a least-squares 
method are briefly discussed. Possible improvements in the 
method are discussed but, as of the present, it- appears im- 
provements would be justified only for those cases where the 
basic assumptions are closely approached and where instru- 
ment accuracy is high. The method is equally applicable to 
other problems which can be expressed by second-order 
differential equations. 

Langley Aeronautical Laboratory, 

National Advisory Committee for Aeronautics, 
Langley Field, Va., December 15, 1950. 


APPENDIX 

FREQUENCY-RESPONSE RELATIONS 

In the body of the report the phase angles aDd amplitude ratio were given for only the simplest case. The complete 
frequency -response relations and transfer functions for the equations involving 5 and its derivatives are now presented. 
If D represents the differential operator djdt, then the steady-state response due to a sinusoidal forcing function can he ob- 
tained by substituting iw for D in the transfer function. The following relations were developed by this procedure: 


Equation 

Transfer Function 

Rhase Angle 

Amplitude Ratio 

a + Kra+KtA^KsAS+Kj 

A a 

Ki+KJ) 

tb tan -1 K\d) 11-1 



Aa 

1 A 3 2 +A 4 co 2 

AS 

d°-+KiD+k 2 

Vctj lilil jy- 2 1 

4 K 2 — 

Kz 


\5 

M (a 2 -^) 2 +atv 


A a 

Kz 




Aa 

Kz 

a -\- K\a + K .2 Aq:= K% Ad 

AS 

d-+k x d+k 2 

^ 5 — tan K.% — c *)2 



AS 

V(A 2 - W 2 ) 2 + A! 2 W 2 









6-^ K 2 Ad=K^Ad-\'K^ ^ Ad dt 

AO 

K,D+K S 

A, ton _1 ^ W I tan -1 


A6 


/ Ae 2 +A 5 V 

A5~ 

D(D*+KiD+K 2 ) 

1 Pd: L-clIi. -j--r | UUl 

d ill CO 

Kq 

AS 


V ^[AtV+d^-w 2 ) 2 ] 







H -f- K\ ft Kz All— K-? Ad -f- "I - K$d 

An 

A 7 +A 8 D+A 8 £ 2 


K S d> 

Ai 

l /A 8 V + (A;-A 5 o /) 2 

A5~ 

D ? + K:D + K 2 

— tan ^_ w2 +tan R _ 

-A „« 2 

A 

5 1 

V (Ki-oty+KW 

7 i -j- K i ib -f- A 2 An = Aj A 5 — K%S 

An 

K : +K S D 


Aged 


An 

1 AsV+Ar 2 

AS 

D’+A.D+A, 

4>JtX k&Tl . rr 2 1 ^111 

0 ii 2 — or 

a 7 


AS 

V (A 2 — w 2 ) 2 -)- Ariw 2 

n +K ] h+ K 2 A n. = A? A 5 

An 

Kj 

, , 1 —Kid) 



An 

1 k 7 

A5 

W+KiD+k, 

^j-^in Ki—tf 



AS 

I V(A 2 -a > 2 ) 2 + AT 2 co 2 
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